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(Dated:) 

We propose a variational approach to study renormalized phonons in momentum conserving nonlinear lattices 
with either symmetric or asymmetric potentials. To investigate the influence of pressure to phonon properties, we 
derive an inequality which provides both the lower and upper bound of the Gibbs free energy as the associated 
variational principle. This inequality is a direct extension to the Gibbs-Bogoliubov inequality. Taking the 
symmetry effect into account, the reference system for the variational approach is chosen to be harmonic with 
an asymmetric quadratic potential which contains variational parameters. We demonstrate the power of this 
approach by applying it to one dimensional nonlinear lattices with a symmetric or asymmetric Fermi-Pasta- 
Ulam type potential. For a system with a symmetric potential and zero pressure, we recover existing results. 

For other systems which beyond the scope of existing theories, including those having the symmetric potential 
and pressure, and those having the asymmetric potential with or without pressure, we also obtain accurate sound 
velocity. 

PACS numbers: 63.20.-e, 63.20.Ry, 45.10.Db, 44.05.+e 


I. INTRODUCTION 

Fleat conduction in low dimensional anharmonic systems 
has attracted considerable interest in recent years [1-3]. 
Phonon, as the predominant heat carrier in insulating materi¬ 
als, undoubtedly lies in the heart of heat conduction. However, 
phonon bears a solid basis only in harmonic systems. The 
intrinsic nonlinearity in anharmonic systems will inevitably 
affect the behavior of phonon. Therefore, understanding the 
phonon properties in anharmonic systems represents an im¬ 
portant question in the heat conduction field. 

A renormalized phonon (r-ph) picture was then put forward 
independently by several groups using varying techniques [4— 
11]. Within the scope of this picture, one can successfully 
interpret and understand a wide range of physical phenom¬ 
ena, including a theoretical description of the sound velocity 
[12] as well as scaling laws of thermal conductivity re(T) with 
temperature T [10, 13]. 

However, the generality of existing, state-of-the-art quasi¬ 
harmonic theories is limited. They were found to provide 
inaccurate predictions on the sound velocity of a nonlinear 
lattice with an asymmetric inter-particle potential [14]. Re¬ 
cently, a numerical method which aims to justify the validity 
of phonon concept in nonlinear lattices was proposed [15]. 
The existence of phonon modes in nonlinear lattices with 
asymmetric potentials is confirmed for a wide range of param¬ 
eters. Hence, a unified quasi-harmonic theory which extend to 
cover the general cases beyond harmonic is desirable. 

In statistical mechanics, variational approaches are often 
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used to obtain approximate information of a nonlinear system 
via a reference system and an associated variational principle 
[16]. The reference system whose properties can be easily ob¬ 
tained contains several variational parameters. In a variational 
scheme, an optimal reference system can be obtained by vary¬ 
ing those parameters such that bounds of the variational prin¬ 
ciple go to a relative minimum or maximum. 

In this paper, we develop a variational approach to study 
phonons in nonlinear lattices. We choose a general harmonic 
system as the reference system for this purpose and regard 
the so-obtained optimal harmonic system as an approximation 
to the original system from which we identify the properties 
of r-phs. We consider applications to one dimensional (ID) 
nonlinear lattices with either symmetric or asymmetric inter¬ 
particle Fermi-Pasta-Ulam (FPU) potentials [17] and demon¬ 
strate the power of our approach by comparing with molecule 
dynamics (MD) results. 

This paper is organized as follows. In Sec.II, we introduce 
our variational approach. In Sec.Ill, we present numerical de¬ 
tails. In Sec.IV and Sec.V, we study nonlinear lattices with 
symmetric and asymmetric FPU potentials, respectively. Fi¬ 
nally, we briefly summarize the work in Sec. VI. 


II. THE VARIATIONAL APPROACH 
A. The variational principles 

Our current investigation aims to develop a quasi-harmonic 
theory applicable for nonlinear lattices with asymmetric po¬ 
tentials. One of the attractive features of such systems is the 
existence of nonzero internal pressure due to the asymmetry of 
the potential. Then from a statistical mechanics point of view, 
it’s convenient to use the language of the isothermal-isobaric 
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ensemble which maintains constant temperature T, constant 
particle number N and constant pressure P (so-called NPT 
ensemble) to describe them [18]. Their equilibrium thermo¬ 
dynamic properties can thus be well determined by the Gibbs 
free energy 


G = 


— (3 t 1 In 


j e ~MH + PL) d(ldp 


(1) 


where H, P and L denote the Hamiltonian, pressure and 
volume (or length in ID cases) of the system, respectively, 
Pt = 1 /T is the inverse temperature (we set fcs = 1), q 
and p are short for the products of all the coordinates and mo¬ 
menta of the system, respectively. The volume L is a function 
of q only. Introducing 


H = H + PL , (2) 

whose ensemble average is just the enthalpy of the system 
[19]. The corresponding probability measure in phase space 
reads 


p(q-p) 


f e~P Tn dqdp 


= e -f>T{H-G) 


(3) 


A variational principle is an inequality satisfied by the phys¬ 
ical quantity in which we may be interested [16]. Hence we 
should look for inequalities for the Gibbs free energy G. To 
do this, we introduce a reference system with a Hamiltonian 
H 0 and prepare the original system and the reference system 
at same pressure. Then following this spirit, we integrate both 
sides of the following equality 

p( q,p) = p 0 (q,p)e- /3T(W - Wo) - fe(G °- G) (4) 


Fq, respectively, the enthalpy goes to the energy, the probabil¬ 
ity measure Eq. (3) should also be replaced by the canonical 
measure 


p c (q,p) 


e~^ H 

f e~P TH dqdp 


= e -P t(H-F). 


( 8 ) 


Then we found that the inequality Eq. (6) recovers the well 
known Gibbs-Bogoliubov (GB) inequality [21] 


F < F 0 + (H - H 0 ) p . , (9) 

where Fq is the Helmholtz free energy of the reference sys¬ 
tem and (-)pc stands for the ensemble average with respect to 
the canonical measure p§(q, p) [= e -PT(H 0 -F 0 )j And the 
inequality Eq. (7) goes to 


F > F 0 + (H — H 0 ) pc (10) 

with (-)pc the ensemble average with respect to the probabil¬ 
ity measure p c { q, p) [c.f. Eq. (8)]. As a lower bound of the 
Helmholtz free energy F [16, 22], it’s little applied since the 
ensemble average in it can not be evaluated analytically. How¬ 
ever, it is more accurate than the upper bound in determining 
the free energy of solids [23, 24]. As can be seen later, in 
our case for determining the sound velocity for anharmonic 
lattices, it is still the lower bound that gives better prediction 
comparing to the upper bound. 


B. The nonlinear systems 

We consider ID momentum conserving nonlinear lattices 
described by the general Hamiltonian [25] 


over the whole phase space to yield 

1 = e ~0T(G o -G) . / ) ( 5 ) 

' ' Po 

where Gq is the Gibbs free energy of the reference system and 
Hq = Hq + PLq, (-) Po denotes the ensemble average under 
the probability measure po(q^P) = 

Taking logarithm over both sides of the above equation and 
using the Jensen’s inequality for exponential functions, i.e., 
(e x ) ^ e^(following [20]), we have 

G < G 0 + (n-H 0 ) po . (6) 

By switching the role of the nonlinear system and the refer¬ 
ence system in Eq. (4), we obtain 


G > G 0 + (H - Ho) P (7) 

with (-)p the ensemble average with respect to the probability 
measure p(q, p) [c.f. Eq. (3)]. These two inequalities [Eqs. 
(6) and (7)] give the upper and lower bound of G, respectively. 

If the pressure vanishes, i.e., P = 0, the Gibbs free ener¬ 
gies G and Gq reduce to the Helmholtz free energies F and 


N 


*= £ 


— Qn-l) 

2m 


( 11 ) 


where N is the particle number, p n denotes the momentum of 
n-th particle, q n = x n nr denotes the displacement of n-th 
particle from its equilibrium position nr with x n the abso¬ 
lute position and r the equilibrium distance for the interaction 
bond, and V represents the inter-particle potential. For brevity 
and without loss of generality, we take m = 1 and r = 1 as 
the unit of mass and length, respectively. 

The average lattice spacing a is then given by 1 + (q n — 
Qn-i)p , and the average lattice length L = (L) p equals Na 
for an A’-particle lattice. We say that the lattice is at its natural 
length if a equals r [=1], namely, L = N for an A r -particlc 
lattice. The average length can be changed to other values by 
applying pressure. 

Furthermore, we introduce S n = q n — q n -i- If the potential 
satisfies V(6 n ) = V(—d n ), we refer it to a symmetric poten¬ 
tial, otherwise we call it an asymmetric one. In terms of 5 n 
and p n , the equations of motion (EOMs) read 

3n — Pn -\-1 Pn ? 

Pn = V'(5 n ) - V'(6 n - 1), 


( 12 ) 

(13) 
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where the dot and the prime denote the time and space deriva¬ 
tive, respectively. 

Specifically, in this work, we focus on the FPU lattices with 
the inter-particle potential [17] 

V(S n ) = \f n + + fa, (14) 

which has become an archetype ID nonlinear system in sta¬ 
tistical mechanics [26, 27]. We call the lattice with a = 0 the 
FPU-/3 lattice. Otherwise, we refer it to the FPU-a/3 lattice. 
For simplicity but without loss of generality, we only study 
the case of /3 = 1 in this paper. It is evident that the FPU-/3 
lattice has a symmetric potential while the FPU-a/3 lattice has 
an asymmetric one. Therefore, an FPU-/3 lattice with natural 
length has zero pressure, while an FPU-a/3 lattice with natu¬ 
ral length has a nonvanishing pressure due to the asymmetry 
of the potential. Nevertheless, our variational principles [Eqs. 
(6) and (7)] enable us to study them within the same theoreti¬ 
cal scheme. 


C. The harmonic reference system H o 


Phonons bear a solid basis only in harmonic systems, in 
order to develop an effective phonon theory for nonlinear lat¬ 
tices, we should choose a quadratic Hq as a reference sys¬ 
tem. Moreover, the present work investigates nonlinear lat¬ 
tices with asymmetric inter-particle potentials. Taking those 
into consideration, we consider a harmonic reference system 
in the following form 


N 


H ° = E 


Y + WO 


Vo(5 n ) = y (S n -d) 2 , 


(15) 

in which K and d are variational parameters with K being 
the effective elastic constant and d quantifying the degree of 
asymmetry of the potentials. Note that by setting d = 0, we 
can use Hq to deal with nonlinear systems with symmetric 
potentials as well. 

The dispersion relation of this reference system is given by 


Wfc 


2\fK sin y, 


(16) 


where k is the wavenumber. It can be regarded as an approxi¬ 
mate dispersion relation for the r-phs in nonlinear lattices [15]. 
The corresponding sound velocity reads 


c s = \[Ka. 


(17) 


The variational method is then to select the optimal refer¬ 
ence systems with parameterized Hamiltonian Hq that either 
minimize the right hand side (r.h.s.) of Eq. (6) or maximize 
the r.h.s. of Eq. (7). This strategy then provides approxima¬ 
tions for the Gibbs free energy and allows us to find optimal 
Hq for the system H. The process is going to be detailed in 
the next subsection. 


D. Determining the optimal harmonic systems 


Note that dqidq -2 ■ ■ ■ dq jv = d5\d52 • • • in the large N 
limit (the Jacobian is unity), then for a system with the proba¬ 
bility measure Eq. (3), we have 

N 

P(q^P) = n P*( S n,Pn), (18) 

71=1 


where p s denotes the marginal phase space distribution for the 
single site variables S n and p n of n-th particle [28, 29] 


Ps(5,p) = ^exp<( Pt 


P- + V (5) + P6 


(19) 


with z the corresponding partition function 


z= exp { —/3 t 


y + V(S) + P8 


d5dp. (20) 


Such a result can be readily tested using MD simulations, see 
Sec. III. 

For a system with sufficiently large particle number N, the 
Gibbs free energy can be expressed as 

G = Ng= -Nfc 1 In z, (21) 

where g = —3 ^ 1 In z is the Gibbs free energy per particle. So 
the Gibbs free energy of the reference system Hq reads 

G 0 = -Np^lnzo (22) 


with zo the partition function determined by Eqs. (15) and 
(20). It can be evaluated analytically, which gives 


27T 


ZQ = 


/3Tv/ff" 


exp 


-h\pd-H 


(23) 


Using the single-site distribution Eq. (19), and noticing that 
the length L and Lq have a same expression, namely, ]T/ n ( 1 + 
S n ), the averages in Eqs. (6) and (7) read 


(H - H 0 ) PO = N ■ {V{6) - V 0 (S)) p o ; (24) 

(H - H 0 ) p = N ■ (V(6) - V 0 (6)) ps , (25) 

where (■) 0 and (•) denote single-site averages with respect 
to the probability measure Eq. (19) with the potential being 
Vo and V, respectively. We will use these symbols in the rest 
of the paper. 


1 . The upper bound harmonic system 


Firstly, we focus on the upper bound of G [Eq. (6)]. Min¬ 
imizing it with respect to the variational parameters K and d, 
i.e.. 


-^[G 0 + (H - Hq) P o ] - 0, 
tt-[Go + (H - Hq) P o \ = 0, 


(26) 

( 27 ) 
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we obtain the following coupled equations: 


du 

K v 


((v-iwos + i)^ 

(y ^o) p0 

_ (V-Vo)# _ 

Pt ({V -V 0 ){S-du + ^) 2 ) p0 


(28) 

(29) 


It can be readily tested that the second order derivatives are 
positive at (du, Ku) so that this set of solution indeed gives 
the minimum of the upper bound. Thus Ku and du determine 
the optimal upper bound of the Gibbs free energy G and the 
corresponding optimal Harmonic reference system (denoted 
as the UH system). 



FIG. 1: Average lattice length L as a function of pressure P. The 
dashed-dotted line indicates that for the ID FPU-/3 lattice with T = 
1. The solid line denotes that for the ID FPU-a/3 lattice with T = 1 
and a = 1. 


2. The lower bound harmonic system 

Now we turn to the lower bound of G [Eq. (7)]. Similarly, 
we maximize it with respect to the variational parameters K 
and ri, the optimal solution reads 

dL = (d)p s + — , (30) 

Kl = fhW) P .-m P .n (31) 

We have also checked that Kl and <7/, indeed give the max¬ 
imum of the lower bound, thus them uniquely determine the 
optimal lower bound of the Gibbs free energy G and also the 
corresponding optimal harmonic reference system (denoted as 
the LH system). Note that Ku does not rely on d / and is com¬ 
pletely determined by Eq. (31). 

The LH system and the UH system can be regarded as ef¬ 
fective descriptions of the original nonlinear system. If we use 
these two optimal systems to construct the effective theory of 
r-phs in the original nonlinear system, then the dispersion re¬ 
lation should follow Eq. (16) with K given by Kl or Ku- 
The accuracy of the results can be evaluated by comparing the 
calculated sound velocity with the predictions using Eq. (17). 
In the following, we present numerical details and then apply 
this strategy to two ID nonlinear lattices, namely, the FPU-/3 
lattice and the FPU-a/3 lattice. 

III. NUMERICAL DETAILS 


average length L equals N(l + (S) Ps ), so a stressless sys¬ 
tem with an asymmetric potential corresponds to an average 
length N(l + (6) p c), where p c s means the canonical single-site 
distribution, i.e., p s [Eq. (19)] with P = 0. Particularly, for 
the FPU-a/3 lattice with zero pressure, the average length is 
smaller than N as can be seen from the figure. 

An NPT system with pressure P can be prepared either by 
applying an external pressure P to the end particles with free 
boundary condition or by fixing the total length L to L(P), 
where L(P ) as a function of P is just depicted in Fig. 1. 
In the present work we adopt the second approach. Specif¬ 
ically, a modified periodic boundary condition is thus used 
to fix the total length at a certain value, namely, qN — q o = 
L — N = N(S) Pe . For such systems, a symplectic integra¬ 
tor SABA 2 with a corrector SABA 2 C [30] is adopted to in¬ 
tegrate the EOMs with a time step h = 0.02 to ensures the 
conservation of the total system energy and momentum up to 
high accuracy 10 -14 in the whole run time, a transient time of 
order 10' is used to equilibrate the system. 


A. Single-site distributions 

We first check the single-site distributions in Eq. (19). After 
the system reaches its thermal equilibrium state, we can obtain 
time series of v n (equals p n ) and 5 n of the n-th particle (n can 
be an arbitrary integer from 1 to N). We then plot their his¬ 
tograms and compare the envelopes against p s = p v ps with 


In this work, we consider systems either with or without 
pressure P. The average lattice length L can be controlled 
by adjusting the pressure, as demonstrated in Fig. 1. It can 
be clearly seen that for the FPU-/3 lattice once the pressure 
is nonzero, the average length of the lattice is away from the 
natural length, i.e., L N. It further indicates that we can 
apply a positive or negative pressure to this lattice system in 
order to compress or elongate its total length. But for sys¬ 
tems with asymmetric potentials, i.e., the FPU-a/3 lattice, the 
pressure would be nonzero at its natural length. Note that the 



where zg = f dS exp[— Pt(\8 2 + § S 3 + \8 A + P<5)]. 

To illustrate the comparison, we simulate the FPU-/3 lattice 
with N = 2048, T = 1, and L = N or L = 1.27V, the 
corresponding pressure is 0 or -0.4309 according to Fig. 1, re¬ 
spectively. A similar simulation is carried out on the FPU-a/3 
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lattice with N = 2048, T = 1, a = 1, and L = 0.8023./V 
or L = N, with the pressure equals 0 or -0.3904, respec¬ 
tively. The results are depicted in Figs. 2 and 3. p v and p$ 
are shown as green and red solid curves in the figures, respec¬ 
tively. Perfect agreements between theoretical curves and the 
MD results are presented, which indicates the validity of the 
single-site distribution in Eq. (19). 


or strong anharmonicity will affect its accuracy. The second 
one is to look at the frequency of the lowest phonon peak of 
the power spectrum of an A r -particlc lattice [14] 

wi = 2— sin^r, (34) 

a N 



FIG. 2: Single-site distributions for the ID FPU-/3 lattice with P = 0 
(upper panel) orP ^ 0 (lower panel). The blue parts in the left 
and right column denote the histogram of <5 and v, respectively. The 
green and red solid lines are p v and pg [c.f., Eqs. (32) and (33)], 
respectively. 



FIG. 3: Single-site distributions for the ID FPU-a/3 lattice with with 
P = 0 (upper panel) or P ^ 0 (lower panel). The blue parts in the 
left and right column denote the histogram of J and v, respectively. 
The green and red solid lines are p v and pg [c.f., Eqs. (32) and (33)], 
respectively. 


B. Sound velocity 


where c s is the sound velocity. This lowest phonon peak can 
always be detected in the power spectrum with high resolution 
(see Fig. 4), thus the sound velocity can be well determined by 
u> i. The third one obtains the sound velocity from the disper¬ 
sion relation [15]. This method, although follows the defini¬ 
tion of the sound velocity, is most computationally expensive 
comparing with the other two, because the calculation of the 
dispersion relation is time consuming. 

Therefore, in this study, we choose the second method to 
obtain MD results of the sound velocity of nonlinear lattices. 
According to the Wiener-Khinchin theorem [32], the power 
spectrum P(ui) can be obtained by doing Fourier transform of 
the velocity autocorrelation of a single particle {v n (t)v n (0)) 
with t = 0, h, 2 h, ■ ■ ■ , t ma x- In order to get smooth phonon 
peaks in the power spectrum, the maximum correlation time 
tmax whose inverse 1 /t max determines the frequency resolu¬ 
tion is set to be 2 M h under the condition that N = 1024. Fig. 
4 presents the power spectrum of the FPU-/3 lattice and the 
FPU-a/3 lattice with a = 1 at their natural length. T = 1 is 
used for both models. It is apparent that phonon peaks with 
the lowest frequencies are well distinguished from the power 
spectrum. The nonlinearity renders its effect in the phonon 
peak broadening. Interestingly, we found that the phonon 
peak of the FPU-a/3 lattice is broader than the one of the 
FPU-/3 lattice, which indicates that the asymmetric potentials 
provide stronger phonon-phonon interaction compared to the 
symmetric ones. 



Now we turn to the calculation of the sound velocity. So far, 
there are mainly three numerical methods which can obtain 
the sound velocity of the nonlinear lattices. The first method 
regards the moving velocity of front peaks of the equilibrium 
spatiotemporal correlation function as the sound velocity [12, 
31]. However, a broadening of front peaks at high temperature 


FIG. 4: Power spectrum P(lo) in the low frequency regime. The 
dashed-dotted line stands for the FPU-/3 lattice at its natural length. 
The solid line represents the FPU-a/3 lattice with a = 1 and natural 
length. For comparison, the lowest harmonic frequency is depicted 
as black solid line. 
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IV. THE FPU-/3 LATTICE 

In the following, we will give a detailed study of the FPU- 
/? lattice by using our variational approach. Being a special 
case of the FPU-a/3 lattice, i.e., a = 0, it has a symmetric 
inter-particle potential. 

A. P = 0 

The existing quasi-harmonic theories of this model all focus 
on this particular case [4, 5, 7-11]. In this part, we not only 
present the comparison between our approach and MD results, 
but also point out the connections between our approach and 
some of the existing quasi-harmonic theories. The first obser¬ 
vation is djj = d,L = 0 [c.f., Eqs. (28) and (30)] by noting 
that the potential of the FPU-/3 lattice is an even function of 
6. Then in the following, we only need to concern about the 
parameter K. 

Firstly, we focus on the UH system. The corresponding 
parameter Kjj [Eq. (29)] can be calculated to yield 

Kl - K v - 3/3^ = 0, (35) 

which has a positive solution 

Ku = 7^(1 + \A + 12/3^ 1 ). (36) 



FIG. 5: Sound velocity as a function of temperature for the ID FPU- 
/3 lattice with zero pressure. The circles are MD results obtained 
from the power spectrum, the solid line the prediction of the LH sys¬ 
tem (\/Kl [Eq. (37)]), and the dashed-dotted line the prediction of 
the UH system {\JKy [Eq. (36)]). The inset presents the relative 
deviation of the sound velocity: a = (c™ — c^ 1D )/c^ ID • 100% with 
c™ and c“ D the theoretical result and the MD result of the sound ve¬ 
locity, respectively. The squares denote c™ = \ /Kl [Eq. (37)] and 
yjKu [Eq. (36)] in red and blue, respectively. 


specific situation, i.e., the lattice length is fixed to be 1.2AT. 
Since L = N( 1 + (S) flg ), we have {S) pg = 0.2. Therefore, 
the pressure P can be obtained by solving 


This result coincides with the prediction of the so called self- 
consistent phonon theory (SCPT) [10] for the FPU-/? lattice, 
since the SCPT based on the GB inequality Eq. (9) as well. 

Then we turn to the LH system. The parameter Kl [Eq. 
(31)] reduces to the simple form 


Kl 


1 


(37) 


The corresponding sound velocity i JKl is exactly the same 
as that predicted by the effective phonon theory (EPT) [8, 12] 
based on the generalized equipartition theorem [33] as well 
as the nonlinear fluctuating hydrodynamics (NFH) using the 
hydrodynamic approximation [29, 34, 35], although the NFH 
is not an effective theory for phonons. 

Results of the sound velocity are illustrated in Fig. 5. The 
excellent agreement between the LH system’s predictions and 
MD results is obvious. While for the UH system, the deviation 
from MD results becomes larger and larger as the temperature 
increases. The reason behind this fact is clear, since high order 
nonlinear terms ignored by the UH system become important 
at high temperature. Hence the LH system can describe r-phs 
in the FPU-/3 lattices without pressure. 


B. P ± 0 

Now we begin to apply our variational approach to the 
FPU-/3 lattice with pressure. The method used to obtain 
nonzero pressure has been discussed in Sec. III. For sim¬ 
plicity, but without loss of generality, here we only study a 


._ f5e-^ v + PS US 

~ fe-fr(v+P6) dS ~ °' 2 ' 


(38) 


The pressure is temperature dependent, which is plotted in the 
inset (b) of Fig. 6. Since the lattice is slightly elongated, the 
pressure should be negative. 

With those values of pressure, we can calculate the aver¬ 
ages in Eq. (31) for Kl numerically and thus obtain the value 
of Kl- As for Kjj, we insert the values of pressure into Eq. 
(28) and Eq. (29), then solve the coupled equations together 
to get the value of Ku- The sound velocity of the LH system 
and the UH system can be obtained by i -JKya and \/ Ky a, re¬ 
spectively. We present results of the sound velocity in Fig. 6. 
It is clearly shows that the LH system still gives better results 
than the UH system, a significant deviation exists between the 
latter and MD results. So the LH system can also describe 
r-phs in the FPU-/3 lattices with pressure. 

The prediction of the sound velocity given by the NFH is 
also shown in Fig. 6. It takes a complex form 


c s /a 


( ±/?r 2 + (V + PS-,V + PS) 

y/K NFH, 


1/2 


(39) 


in which all the (•) is short for {-) Pe , i.e., averages under 
the single site probability measure Eq. (19), and (A; B) = 
(AB) — (A) ( B ). Noticing the EPT and the SCPT are unable 
to deal with systems with pressure, they become invalid in this 
case. 
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FIG. 6: Sound velocity as a function of temperature for the ID FPU- 
/? lattice with nonzero pressure. The circles are MD results obtained 
from the power spectrum, the solid line the prediction of LH sys¬ 
tem WKl [Eq. (31)]), the dashed line the prediction of the NFH 
WKnfh [Eq. (39)]), and the dashed-dotted line the prediction of 
the UH system (y/Ku [Eq. (29)]). Inset (a) presents the relative de¬ 
viation of the sound velocity: a = (c™ — c^ ID )/c^ 1D • 100% with 
c™ and c^ D the theoretical result and the MD result of the sound 
velocity, respectively. The squares denote c™ = \/Kl [Eq. (31)], 
\/Kjj [Eq. (31)], and \[K nfh [Eq. (39)] in red, blue, and cyan, 
respectively. Inset (b) shows pressure P as a function of temperature 
T for ID FPU-/? lattice with L = 1.2 N. 


V. THE FPU-a/3 LATTICE 


When a cubic term is added to the FPU-/3 potential, we get 
the asymmetric FPU-a/? potential [Eq. (14)]. The potential 
with different a is plotted in Fig. 7. It is clearly seen that 
the potential becomes more asymmetric as a increases. When 
a exceeds 2, dV/dS = 0 begins to have two solutions, the 
potential tends to become a double-well one which is out of 
the scope of the present investigation, we can see that from 
the form of the potential with a = 2.1 in the figure. The 



A. P = 0 


We first apply the UH system to the FPU-a/3 lattice. The 
corresponding parameters du and Kjj now satisfy 


2 dfr + dfr 


4 

1 


3\Kl 


1 

/3t 


( djj -p cxd^ 


PtKu 
(a + 3 du) = 0, 

(1 + 2adjj 4 


PtKu 


du) 


3 dfj) — 


3\Kl 


(40) 


= 0.(41) 


We note that du is no longer zero due to the asymmetric po¬ 
tential V. The two equations are coupled, we should solve 
them together to get the value of Ku- 

We then utilize the LH system to investigate the lattice. The 
parameters of the LH system are given by 


d-L = (d) pi , 

(42) 

Kl ~ pT[(s 2 ) P o~((d) p oyy 

(43) 


Similarly, dz, is nonzero for an asymmetric potential. Note 
that the EPT still predicts the effective force constant as [8] 


Kept 


1 

JiWbi 


(44) 


with V now the FPU-ci/3 potential Eq. (14). It is evident that 
the denominator of Kl and Kept take the form of the vari¬ 
ance and the second moment of <5, respectively. This discrep¬ 
ancy results from the asymmetry of the potential, since (6) p c 
vanishes for symmetric potentials. As can be seen later, such 
a correction improves the accuracy of the results significantly. 

We consider lattices at a fixed temperature T = 0.5 and 
vary a from 0.2 to 2. Results of the sound velocity are illus¬ 
trated in Fig. 8. The deviation of predictions of the EPT from 
the measured value is clearly revealed as a increasing. This 
demonstrates the invalidity of the EPT for the strong asym¬ 
metric cases. While the LH system gives accurate results no 
matter how large the asymmetry is. It thus shows the signif¬ 
icance of the correction term in the denominator of Kl [Eq. 
(43)]. Still, we observe that the LH system gives much better 
results than the UH system. So we can regard the LH system 
as an effective description of r-phs in the FPU-a/3 lattice with¬ 
out pressure. Meanwhile, an agreement between the NFH and 
the LH system renders a fact that hydrodynamic approxima¬ 
tions can also be applied to long-wavelength r-phs in lattices 
with asymmetric potentials. 


B. P^0 


FIG. 7: Potential V(S) of the ID FPU-a/3 lattice with varying a. 


inherent asymmetry can induce thermal expansion, which has 
been shown to play a significant role in heat conduction [36- 
39]. In the following, we will study this lattice by using our 
variational approach. 


Now we consider the FPU-a/3 lattice with pressure. For 
simplicity. We only investigate the system at its natural length. 
The value of the pressure can be obtained by solving 


fie-M v + p Qd6 

Wp« “ J e -PT(v+P6)dS ~ U ‘ 


( 45 ) 
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FIG. 8: Sound velocity as a function of a for the ID FPU-a/3 
lattice with T = 0.5 and zero pressure. The circles are MD re¬ 
sults obtained from the power spectrum, the dashed line the predic¬ 
tion of the NFH (\/Knfh with P = 0 [Eq. (39)]), the dashed- 
dotted line the prediction of the UH system (y/Ku [Eq. (41)]), the 
red solid line the prediction of the LH system (\[Kl [Eq. (43)]), 
and the maroon solid line the prediction of the EPT (n/Kept [Eq. 
(44)]). The inset presents the relative deviation of the sound velocity: 
a = (c™ — cT)/cr ' 100% with c™ and c“ D the theoretical result 
and the MD result of the sound velocity, respectively. The squares 
denote c™ = \JKl [Eq. (43)], \JKp [Eq. (41)], \JKnfh with 
P = 0 [Eq. (39)], and \] Kept [Eq. (44)] in red, blue, cyan, and 
maroon, respectively. 


Results of pressure as a function of temperature are plotted 
in the inset (b) of Fig. 9. The absolute value of pressure 
increases as the temperature increases because of the inter¬ 
particle interaction is stronger. At the low temperature regime 
where thermal excitations dwell the very bottom of the po¬ 
tential, particles can not feel the asymmetry of the potential 
strongly, so the pressure approaches zero as the temperature 
tends to zero. 

With these values, Kl can be easily calculated from Eq. 
(31). The value of Kjj can be obtained by solving Eqs. (28) 
and (29) together. Since the existing quasi-harmonic theo¬ 
ries fail in this model, only our results and the NFH’s will be 
shown. 

The result for the case a = 1 is illustrated in Fig. 9. From 
this figure, we see that the LH system works better than the 
UH system as usual. 

We further investigate the FPU-a/3 lattices with varying a 
by using the LH system. The temperature T is fixed to be 0.5. 
From the Fig. 10, it is apparent that the discrepancy between 
the theoretical prediction and MD results tends to increase as 
a increases. One possible reason is that the Gibbs free energy 
given by the lower bound depart from the exact result signif¬ 
icantly at large a (we can see that from the inset in Fig. 10). 
This fact indicates that the first cumulant approximation we 
adopted in deriving the inequalities is not enough. Higher or¬ 
der contributions to the Gibbs free energy must be considered 
in order to improve this method [40, 41]. 



FIG. 9: Sound velocity as a function of temperature for the ID FPU- 
a/3 lattice with a = 1 and nonzero pressure. The circles are MD re¬ 
sults obtained from the power spectrum, the solid line the prediction 
of LH system (s/Kl [Eq. (31)]), the dashed line the prediction of the 
NFH (\/Knfh [Eq. (39)]), and the dashed-dotted line the prediction 
of the UH system (y/Ku [Eq. (29)]). Inset (a) presents the relative 
deviation of the sound velocity: a = (c™ — c^ D )/c^ D • 100% with 
c™ and the theoretical result and the MD result of the sound 
velocity, respectively. The squares denote c™ = \/Kl [Eq. (31)], 
y/Kl 7 [Eq. (31)], and sfK nfh [Eq. (39)] in red, blue, and cyan, 
respectively. Inset (b) shows pressure P as a function of the temper¬ 
ature T for the FPU-a/3 lattice with a = 1 at its natural length. 



FIG. 10: Sound velocity as a function of a for the ID FPU-a/3 lat¬ 
tice with T = 0.5 and nonzero pressure. The circles are MD results 
obtained from the power spectrum, the solid line the prediction of 
V Kl [Eq. (31)]. The inset presents the Gibbs free energy per par¬ 
ticle. The dashed line is the exact Gibbs free energy g according to 
Eq. (21), the solid line is the prediction of the lower bound according 
to Eq. (7). 

VI. SUMMARY 

In summary, we have presented a variational approach to 
study renormalized phonons in momentum conserving non¬ 
linear lattices. Specifically, we obtain two optimal harmonic 
reference systems, namely, the LH system and the UH sys¬ 
tem, by optimizing the two bounds of the Gibbs free energy 
of the nonlinear system. These optimal systems which deter¬ 
mine the properties of renormalized phonons can be regarded 
as optimal harmonic approximations of the nonlinear system. 
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This method has been applied to lattices with either sym¬ 
metric or asymmetric potentials, with or without pressure. For 
lattices with symmetric potentials, such as the FPU-/? lattice, 
our variational approach gives two optimal and symmetrical 
reference harmonic lattices. In the zero pressure case, it is 
very interesting to find that the sound velocity derived from 
the LH system reproduce the previous theoretical results from 
the effective phonon theory and nonlinear fluemating hydro- 
dynamic theory, and the sound velocity derived from the UH 
system recover the previous theoretical results from a self- 
consistent approach. 

For lattices with asymmetric potentials, such as the FPU- 
a/3 lattice where the existing effective phonon theory fails to 
predict the sound velocity, our variational approach can also 
give accurate predictions. In particular, our approach reveals 
the reason why the effective phonon theory cannot predict the 
correct sound velocity. 

In all the cases, the LH system works better than the UH 
system, since the ensemble averages in the former are eval¬ 
uated under the probability measure of the nonlinear system, 
it simultaneously takes nonlinear contributions into account 
compared with the latter. Therefore, the LH system can be 
treated as an effective theory of renormalized phonons in var¬ 
ious momentum conserving nonlinear lattices. A deviation 
between the prediction of our variational approach and the 
MD results has also been found for stronger asymmetry and 
nonzero pressure situation where the underlying mechanism 
needs further investigations. 

Compared with the existing quasi-harmonic theories, this 
approach can be used to investigate nonlinear systems with 


asymmetric potentials. We owe this ability to two aspects. 
One is the choice of an parameterized asymmetric harmonic 
potential with a parameter d which can quantifies the degree 
of asymmetry of the potentials. The other is the consideration 
of the Gibbs free energy instead of the Helmholtz free energy, 
which enables us to deal with systems with pressure. 

We also found that the NFH gives satisfactory predictions 
of sound velocity in all those cases. The theory focuses on 
a hydrodynamic description of nonlinear lattices, not a con¬ 
struction of an effective theory of the renormalized phonons 
in nonlinear lattices. Those agreements means that the hydro- 
dynamic approximation really captures essential features of 
the systems in the long-wavelength limit. However, for sys¬ 
tems without momentum conservation, the concept of sound 
velocity is invalid, the NFH may not be able to give the in¬ 
formation of the renormalized phonons, although it can still 
describe correlation functions of those systems [42]. While 
our variational approach can be uesd to investigate phonon 
band gap and phonon dispersion of those systems [43]. 
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